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We derive the normal modes for a rotating Coulomb ion crystal in a Penning trap, quantize the 
motional degrees of freedom, and illustrate how they can by driven by a spin-dependent optical 
dipole force to create a quantum spin simulator on a triangular lattice with hundreds of spins. 
The analysis for the axial modes (oscillations perpendicular to the two-dimensional crystal plane) 
follow a standard normal-mode analysis, while the remaining planar modes are more complicated 
to analyze because they have velocity-dependent forces in the rotating frame. After quantizing the 
normal modes into phonons, we illustrate some of the different spin-spin interactions that can be 
generated by entangling the motional degrees of freedom with the spin degrees of freedom via a spin- 
dependent optical dipole force. In addition to the well-known power-law dependence of the spin-spin 
interactions when driving the axial modes blue of phonon band, we notice certain parameter regimes 
in which the level of frustration between the spins can be engineered by driving the axial or planar 
phonon modes at different energies. These systems may allow for the analog simulation of quantum 
spin glasses with large numbers of spins. 

PACS numbers: 37.10.Ty , 03.67.Lx, 75.10.Jm, 75.10.Nr 



I. INTRODUCTION 

Richard Feynman motivated the idea of using analog 
quantum simulation as a means to describe the behavior 
of complex quantum- mechanical systems [1 . This idea 
has been experimentally realized for the transverse-field 
Ising model in cold ion traps in one-dimension [IHS], in 
neutral atoms driven to a nonequilibrium state [7], and 
with preliminary results also available for cold ion traps 
in two-dimensions [8] . The tunability and control of cold 
atom and cold ion systems are unmatched in traditional 
(real-material) condensed matter systems and therefore 
they are an ideal platform for examining idealized con- 
densed matter systems and searching for nontrivial quan- 
tum ground states. In addition, the isolation of the sys- 
tem from the environment with long decoherence times 
(in milliseconds) also offers tremendous opportunities for 
the study of non-equilibrium coherent many-body dy- 
namics in real time such as thermalization [9j [10] and 
quench dynamics [TTl [12] in cold atoms. 

Cold ions in traps have been proposed theoretically as 
potential emulators for effective spin models [13] . For the 
case of linear ion chains in Paul traps, the experimental 
protocol has been well established [IHS]. However, there 
are still difficulties in scaling up Paul-trap based systems 
to more than approximately sixteen ions while main- 
taining good quantum control, although next generation 
traps are planned for up to potentially 100 ions. Theoret- 
ically, phonon modes in linear Paul traps have been well 
understood for a decade [14] and are essential for the real- 
ization of the spin models in the system (although effects 
of micromotion on two-dimensional crystals in Paul traps 
are more complicated [E]). On the other hand, ion crys- 
tals in Penning traps can be easily stabilized with several 



hundreds of ions in a two dimensional structure and can 
potentially simulate quantum phases beyond what mod- 
ern computers can simulate. Recently, a feasibility study 
has shown that it is possible to generate spin-spin inter- 
actions in these systems that decay like a tunable power 
law for long distances [8 . Similar to cold ions in linear 
Paul traps, a theoretical understanding of phonon modes 
is needed as a prerequisite for the adiabatic state evolu- 
tion of these systems. 

In this paper, we provide a complete theory for the 
oscillatory normal modes in a Penning trap, tuned to pa- 
rameters similar to those used in current experiments. 
We then quantize these phonons and show how they can 
be used with a spin-dependent optical dipole force to gen- 
erate effective spin- spin interactions. Recently, others 
have shown how to find equilibrium positions for similar 
trapping potentials [16 , and have evaluated the phonon 
spectra using different methods from the ones we em- 
ploy HZ]. 

The organization of the paper is as follows. In Sec. 
II, we formulate the theory for the normal modes of cold 
ions in Penning traps. We determine the ground state of 
the crystal configuration by finding the static equilibrium 
of the system in the rotating frame of a rotating crys- 
tal. The normal modes are then found by harmonic ex- 
pansion around the equilibrium state. Because the Pen- 
ning trap creates a rotating crystal, the normal modes 
must be solved in the rotating frame, and the longitu- 
dinal (planar) modes require the solution of a quadratic 
eigenvalue problem because they have velocity-dependent 
forces. The phonon Hamiltonian is then derived by quan- 
tizing the normal modes. In Sec. II, we also discuss 
quantum simulation based on spin-dependent forces aris- 
ing from laser-ion dipole interactions. The effective Ising 
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models mediated by the axial phonon modes and the pla- 
nar phonon modes are also described. In Sec. Ill, we 
show the numerical results for phonon frequencies and 
for spin-spin interactions for parameters relevant for cur- 
rent experiments. In Sec. IV, we provide our conclusions. 




FIG. 1: Electrostatic potentials provided by electrodes in 
a Penning trap have contributions from both the end-cap 
electrodes [with the trapping potential (/)£;(r) = Vqz'^ that 
pushes the atoms towards z = 0] and the cylindrical electrodes 
[with the radial quadratic potential (/)c(r) = — ^(x^ +2/^)5 
which tends to push the ions out radially]. In addition, 
the repulsive Coulomb potential between the ions tends to 
destabilize the system in the trap. A static magnetic field 
B = Bzz{Bz > 0) is thus applied along the axial direction 
to provide the radial confinement of the ions. To lock the 
rotational angular frequency of the ions at a specific rota- 
tional speed, a time-dependent clockwise quadrupole poten- 
tial (l)w{t) = Vwr^ cos2{0 -\-Qt) (Q > 0) is applied to the ions 
through ring electrodes so that the steady state of the ions 
(with a rigid body rotational speed Q) can be phase locked. 
Note that the rotating quadrupole potential is well imple- 
mented by ring electrodes with just six sectors (as employed 
in the NIST Penning traps (81 [TBI II9] ) . 



of ^Be^ ions, localized in a two-dimensional plane. The 
setup of cold ions in Penning traps is briefly illustrated 
in Fig. 1. Further details can be found in the published 
literature [8, 18 . In current experiments, there also is a 
small amount of impurities such as BeH^ which forms via 
collision of the Beryllium ions with hydrogen molecules 
that are in the background gas. Hence a real system will 
have different mass particles (which will tend to separate 
from one another as the centrifugal force pushes the heav- 
ier particles to the outer regions of the crystal). In this 
paper, we discuss the clean limit only, and ignore such 
defects. The effect of defects will be studied elsewhere. 

Given the above considerations, we form the ion La- 
grangian in the laboratory frame of reference as 
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^m\rj\ 



(1) 



in which N is the total number of ions, e is the positive 
unit charge of and m is the mass of a ^Be^ ion, = 
(pj^Oj^Zj) is the ion spatial configuration in cylindrical 
coordinates, ^j(t) is the total scalar potential acting on 
the jth ion and Aj = ^(B x r^) is the vector potential 
in the symmetric gauge for the axial magnetic field B = 
BzZ {Bz > 0). The potential includes the harmonic 
trapping from the electrodes, the rotating wall potential, 
and the Coulomb interaction between the ions. It satisfies 
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(2) 

in which Vb is the amplitude of the static potential from 
the Penning trap electrodes, Vw is the amplitude of the 
rotating wall potential, 1] > is the rotating wall an- 
gular frequency ft = Qz^ is the Coulomb force con- 
stant, and Tkj = |r/e — rj\ is the inter-particle distance 
between the A:th and jth ions. With the application of 
the quadrupole rotating wall potential ^vF(t), the ion po- 
tential <j)j{t) is time dependent in the laboratory frame 
of reference. However, in the co-rotating frame with the 
angular speed 1], the corresponding effective potential 
becomes time independent and the problem can be un- 
derstood as an equilibrium problem [23 . The geometric 
configuration of the ions under rotation in the steady 
state is simplified by finding the (static) equilibrium ion 
configuration in the rotating frame. In general, the co- 
ordinate transformation between the rotating frame and 
the laboratory frame is described by the 0(2) rotational 
operator Rz{0{t)) with the angle of rotation 6{t): 
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II. THEORETICAL FORMULATION 

We consider hundreds of spins, realized via the two 
hyperfine states pSi/27^j = 1/2), pSi/27^j = 



where the superscript R denotes the rotating frame of 
reference. The operator Rz{0{t)) satisfies 



(4) 
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with 0{t) = —Qt. Accordingly, in the rotating frame, the 
rotating wah potential does not depend on time and 

- Ei Vw{xf 



is given by 



y 



by direct substitu- 



J 



tion. As a consequence, the effective potential energy in 
the rotating frame is then given by the expression 
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which now includes the potential terms from the centrifu- 
gal and Lorentz forces as well. 

Ions trapped in a Penning trap do not always crystal- 
lize into a two-dimensional plane. The following approx- 
imate confinement criterion has to be satisfied to guar- 
antee it is energetically favorable for ions to stay in the 
plane with z = 0: 



26^0 



> 1. 



(6) 



In addition, when the rotating wall potential vanishes 
{Vw 0), the system is confined in the plane only if 
the radial potential is attractive. In this case, the con- 
finement criterion for zero rotating wall potential is given 

by 



(32 = eB^n - m9f - eVo > 0. 



(7) 



For nonzero rotating wall potentials, the ion saddle po- 
tential caused by the rotating wall potential Vw breaks 
rotational symmetry and leads to a confinement of ions 
only when the coefficients of the trapping potential along 
the weak trapping axis is positive. In other words, the 
following criterion has to be satisfied to have a confined 
equilibrium state for all ions: 



■ mn^ - eVo) - e\Vw\ > 0. (8) 



The above heuristic criteria narrow down the relevant 
parameter regime in which a stable planar ion crystal 
formation is feasible for quantum simulation. Of course, 
the accurate quantitative prediction of the stability of 
the ion crystals must also include the Coulomb interac- 
tion between the ions. The precise criterion can only be 
found by solving for the equilibrium positions and show- 
ing that all normal modes have real frequencies, so that 
the equilibrium is globally stable [20 . 

Based on the invariance of the Lagrangian, we derive 
the Lagrangian in the rotating coordinates from Eq. ([T]) 
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(9) 

in which the second term will produce the Lorentz force 
term with the effective mae; netic field B^^f^lO^z for the 



jth ion in the rotating coordinates with the magnitude 
B^^^l^t] = Bz — 2rtm/e which depends on the rotational 
angular frequency 1] > 0. The modification of the mag- 
netic field is due to the Coriolis force for the moving ions 
in the rotating frame. There is no Coriolis force present 
for ions when the ions are in equilibrium in the rotating 
frame but the Coriolis force can have effects on the oscil- 
lating normal mode motion away from equilibrium. The 
centrifugal force originates from the term —^mQ?rf in 
the effective potential energy e(j)f in Eq. dsl. 



A. Stable configurations and normal modes in the 
rotating frame 

It has been well established experimentally that ions 
in Penning traps reach different static equilibrium states 
in the rotating frame of reference (steady state in the 
laboratory frame of reference) for different values of the 
adjustable rotating angular frequency given by a rotating 
wall potential (j)w [22] • Therefore, we can find the stable 
spatial configuration of the ions by determining the local 
minima of the effective potential energy in the rotating 
frame of reference. In general, such a minimization prob- 
lem is difficult to solve for the absolute minimum in two 
or higher dimensions, because many different configura- 
tions can be competitive for the lowest energy state and 
there can be large potential barriers between them. To 
find an experimentally viable solution, we are guided by 
the experiments themselves, which typically observe the 
ions arranged in an approximate triangular lattice in a 
single plane. 

Our strategy is to construct the trial configurations 
based on a "closed shell" construction analogous to find- 
ing stable electronic configurations in an atom as detailed 
in the numerical discussions. Even though we cannot 
guarantee that the stable configuration we find is the 
global minimum, the NIST experimental systems seem 
to find the same local minimum which validates our ap- 
proach. In further support of the strategy, we have suc- 
cessfully predicted the phonon mode spectra and spin- 
spin interaction observed in the NIST experiments [81I2T]. 

In the following, we discuss the collective phonon exci- 
tations about the previously determined equilibrium con- 
figuration. Based on the equilibrium configuration, nor- 
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mal mode dynamics of ions near the equilibrium structure 
can be fully captured by the full system Lagrangian L ex- 
panded around this minimal ion spatial configuration to 
quadratic orders by a Taylor series. We use the follow- 
ing notation for the ion coordinates r j (t) = + SHj (t) 

and for the ion velocities rj{t) = S'Rj{t) in the rotating 



frame. Since the equilibrium configuration is the local 
minimum of the classical action of the system S = J dtL, 
the first order terms in ^Rj(t) and ^Rj(t) do not con- 
tribute to the expansion by construction. Therefore, the 
system Lagrangian can be expanded to quadratic order 
as 



j,k=l L J 1^ -J 



d 
dRk 



(10) 



where the first term Lq is the Lagrangian from the equi- 
librium state and the second term is the Lagrangian Lph 
due to fluctuations away from the equilibrium state. One 
can isolate the planar motions of the ions in the xy plane 
from the axial motions along the z axis based on the 
consideration that the Lorentz force due to the external 
magnetic fields acts only in the crystal plane and the 
potentials are separable in cylindrical coordinates, hence 
there is no harmonic coupling between planar and axial 
degrees of freedom (anharmonic effects can couple nor- 
mal modes together but are ignored here). Therefore, we 
can study the axial and planar modes independently. 

The Lagrangian which governs the axial motions of the 
ions is described as 

N N 

^P/. = 2 E ^^Rj^Rj - 2 E ^jk^Rj^^l (11) 

j=l j,k=l 

in which the symmetric stiffness matrix K^^ is a real 
Hermit ian matrix 
Kj^ are given by 



Hermitian matrix K^^ = K^j and the matrix elements 



Kjk 



where R^j^ 
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|R^— R^ I is the distance between two ions in 



equilibrium (in the rotating frame). Similarly, we derive 
the Lagrangian for the planar normal modes 

N N 

^;. = jEH^i^.i'-jE E Kjk'sR?^R'k 



^ 2 

j=i 



j,k=l a^e(x,y) 



(13) 

in which B^-^-^lTt] = Bz — 2Qm/e is the effective magnetic 
field in the rotating frame with angular frequency Q. The 
real Hermitian stiffness matrix K"^ satisfies the relation 



K^^. We derive the stiffness 



jk kj jk 

matrix K^^ for the planar modes as: 



' m^^ + enB""^ ^ [n] - eVo + 2eVw 



J^jk 



1=1 
2 



DO 5 



kpC 



R%' - 3(xg - 



' mn^ + eVLB'^^^l^] - eVo - 2eVw 



Kfk 



DO 5 



kpC 



N 

1=1 



^jk 



1=1 



DO 5 



j = k,lj^j 



j = kj^j 



(xO-xO)(yO-yO) 



(14) 

where the equilibrium configuration is represented by 
the ion coordinates R^ = {x^^yj^Zj) for j = 1,2,..., A/". 
Notice that the off-diagonal matrix elements K^^ have 
nonzero values, indicating the collective nature of the pla- 
nar motional degrees of freedom, which couple the motion 
in the two coordinate directions together. 



B. Axial phonon modes 

Phonon modes are the quantized modes corresponding 
to the classical normal modes discussed earlier. In this 
section, we solve the classical normal mode problem and 
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then find the quantized Hamiltonian of the axial phonons, 
which is most relevant for the description of conventional 
quantum simulation in a Penning trap. 

By minimizing the action 6 J dtL^^ = 0, the classical 
equation of motion for the ion displacement along the 
axial {-\-z) axis can be written as 



becomes 



N 



m5R] + J2 K^kSRl = 0, i = 1, 2, . . . TV. 



(15) 



k=l 



The axial normal modes are found by direct substitu- 
tion into Eq. (15) of the eigenvector solution SRj{t) = 
bj^ cos[uJzu{t — ^oj] after initial time to where the eigen- 
vectors bj^ are A^-tuples of real numbers with unit norm. 
The following eigenvalue equation must therefore be sat- 
isfied 



N 

E 

fc=i 



[mujlJjk - K^^] hi"" = 0, j> = 1, 2, . . . AT (16) 



where the N axial eigenvectors are the eigenvectors of the 
stiffness matrix K^^ (whose eigenvalues are \zv) but with 
different eigenvalues ujzv = \/ ^zv/'^ (since all eigenval- 
ues \zv are real, the frequencies are either real or pure 
imaginary). Note that this eigenvalue problem is the sim- 
plest quadratic eigenvalue problem, but because one can 
immediately solve it by this mapping onto a linear Her- 
mitian eigenvalue problem, we need not discuss this point 
further here. For the planar phonons, the analysis is more 
complicated, as shown below. 

The planar crystal system is stable against a one-to- 
two plane transition when all the eigenfrequencies uJzv 
are real, which is analogous to case of stability against 
the well-known zigzag transition for ions in linear Paul 
traps. If the eigenfrequencies uOzv are imaginary, then the 
two-dimensional planar equilibrium that was found is an 
unstable equilibrium, and will not form in experiment. A 
detailed numerical analysis of this instability can tell us 
the parameter regimes where a single sheet of ion crystals 
is energetically stable [20 . 

The procedure for the quantization of axial nor- 
mal modes is the standard procedure. One identifies 
the generalized coordinates Qj, and momentum Pj, for 
each phonon mode v as canonically conjugate variables, 
which satisfy the relation given by the Poisson bracket 
{Qvi Pu'} = ^vu'- The quantization of phonon modes can 
be implemented by promoting the relation {Qv^Pv'} = 
5yy' to the commutation relations [Qv^Py'] = ifidyy' for 
operators. To identify the canonically conjugate vari- 
ables for the axial motional degrees of freedom, the La- 
grangian for the axial phonon modes can be recast into 
the form in Eq. ([T7| (as long as the system is in the 
parameter regime in which the planar crystal configura- 
tion is stable with Im[cc;2^] = for all modes v). We let 
SRj{t) = ^j^£,u{t)bj^ be the displacements, where the 
quantum dynamics due to the phonon mode u is given 
by the real generalized coordinate ^u{t)- The Lagrangian 



(17) 



The corresponding generalized momentum for the 
phonon mode u is then given by the expression 



(18) 



As a consequence, the Hamiltonian for the axial phonon 



modes HA is represented by the summation over N in- 



dependent harmonic modes 

N 



A 
ph 



E 



2m 



1 



(19) 



Therefore, the second quantized form of the Hamiltonian 

^ph 



HA for the axial phonon modes simply becomes 



N 



^Ph = ^^^zv{nzv + -), 



(20) 



in which hzv = GL\^azv is the number operator for the 
phonon mode v and the phonon creation and annihilation 
operators are related to the generalized coordinates via 



2% 



2h 



mujzv 
i 

muJzv 



pA 



(21) 



The phonon creation and annihilation operators are re- 
lated to the quantum fiuctuations of ions along the axial 
direction via 



(22) 



It is obvious that the Hamiltonian H^j^ is invariant in 
the rotating frame and the laboratory frame since each 
phonon mode oscillates along the rotation axis z and 
hence they are not infiuenced by the overall rotation of 
the coordinate system. 



C. Planar phonon modes 

So far, the planar phonon modes have not been uti- 
lized for quantum simulation due to the complexity of 
the normal modes as well as the complication introduced 
by the rotation of the ion crystals as observed in the 
laboratory frame. We plan to give a detailed explana- 
tion of the nature of the planar phonon modes and hope 
this will stimulate further insights toward quantum sim- 
ulation. Our procedure is based on previous work on 
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Coulomb crystals in a magnetic field, but not in a ro- 
tating frame [24H27] . Let us expand the Lagrangian 
in terms of the 2N eigenmodes of the hermitian stiffness 
matrix K in Eq. (14), which are labeled by the index 



v with 1 < v < 2N (note that these eigenvectors are 
not the normal mode eigenvectors for the planar motion, 
they are a convenient basis to use for the normal mode 
equations of the planar motion, except for the one case 
noted below). These eigenvectors b^^ satisfy 



E 

k(3 



(23) 



where the spatial indices j = 1, • • • A^, A: = 1, • • • TV are 
the site indices for each ion and a, /3 = x or t/, de- 
note the directions of the spatial displacements of the 
phonons in the plane, mcjQ^ are the real eigenvalues of 
the real-symmetric matrix K for the stable ion configura- 
tion. These eigenvectors do become the eigenvectors for 
the planar modes only in the case when B^f / = (which 
occurs when the rotating frequency is equal to half the 
cyclotron frequency), where the eigenvalue problem sim- 
plifies to the traditional form for a phonon. This fact 
has been used for an analysis of the phonons in a Pen- 
ning trap [17] . The basis vectors b^^ can be chosen to be 
real and they are also orthonormal ^i^^b^^b^^ — ^y^i . 
Therefore, the component of the displacement of the jth 
ion in the ath direction, (t) , can be expanded in this 
basis as 



2N 



SRfit) = J2UtW, j = l,2, 



■N. 



(24) 



where Cu{t) is the dynamic variable. The Lagrangian 
for the planar normal modes in this basis can be rewritten 
following Eq. (13) as 



>2 ^ ,1.2 > 2 



1 



(25) 



where the Einstein summation convention is used for 
repeated indices and the matrix elements Tj^^/ are re- 
lated to the antisymmetric matrix T^^ by the unitary 



transformation 



bTTTu^bi""' in which Tl 



n/3a 



jk 



,q;/3 _ 
jk 

1 and 



^j/c-a/ij — -j-kj (where Cxy — ^yx 
all others vanish) due to the velocity dependent term in 



Eq. (|13|). Based on the Lagrangian L^^^ we can identify 



that THe canonical momentum for mode u is related 
to the mechanical momentum 11^ = mC,v by 



pv 



dL, 



ph 



n'^ - -C'^r T-k'K =11'^- ^C'T.'. 

(26) 

Finally, we arrive at the Hamiltonian for the planar 
modes as 



Notice that the Hamiltonian H^^ is not diagonal in the 
canonical conjugate variables and but is diagonal in 
the variables 11^ and (j^. We need to perform a canonical 
transformation in order to diagonalize the Hamiltonian 
H^^ in the new canonical conjugate variables and we then 
can proceed with the standard quantization procedure 
(as in the axial phonon case). 

Our strategy is to first quantize the Hamiltonian with 
respect to the canonical variables and then choose a non- 
standard linear combination of coordinate and momen- 
tum to create the appropriate raising and lowering oper- 
ators (this strategy was described by Baiko [27 ). Hence 
we use the following fundamental commutation relations 



[C.,C.'] = 0, ] = 0, [CP-^ ] = ihs,. 



(28) 



which lead to the following commutation relations for the 
mechanical momenta 



1 



[H^H^'] = [p^-c.r^f^^] 



,1 



(29) 



in which the antisymmetric property Tj^- 



To 



k -Lf^j u. - -T^/^ is used, 
diagonize the Hamiltonian in the conventional second 
quantized form H^^ = X1a=i ^^a(<^a^a + \) with 
2N collective mode frequencies uo\ > {) for the stable 
ion configurations, we look for the phonon creation 
operator a\ in the presence of the magnetic field as the 
unconventional superposition of the operators H^ and C,y 

a\ = a^lV" + plCv. ax = ^H^ + PTG. (30) 

with a and /3 (possibly complex) numbers. The oper- 
ators ax and must satisfy the commutation relation 
[aA,a]^/] = ^AA' As a consequence, the Hamiltonian H^^ 
is required to satisfy the commutation relation 



ph^ ^a] 



hujxd 



(31) 



with LUx the normal mode frequency for the planar 
phonon. By direct substitution of Eq. (30) into Eq. (31), 



the coefficients and must satisfy the following cou- 
pled eigenvalue equations 



^xa^x = f^x- -T. 

m m 

uJxPl = imujfa^x' 



(32) 
(33) 



Eq. (33) can be solved for /3 in terms of a via 



im^^a\ 



Substituting into Eq. (32), we find the a co- 



efficients satisfy a quadratic eigenvalue problem (QEP) 



- i^xTu 



mujn 







(34) 



H^h = P^Cu - = ^n-^^ + \mu:f(:,\ (27) 



for the eigenvalue uox which is chosen to be a nonnega- 
tive frequency. Since QEP are not so common in physics. 
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we discuss how to solve them in the appendix (see also 
Ref. l28l) . In the simple form in Eq. (34), we can map 



the QEP onto a conventional linear hermitian eigenvalue 
problem in twice as many dimensions, which then allows 
us to use orthonogonality of eigenvectors and complete- 
ness to derive a number of nontrivial identities of the a 
eigenfunctions. Details can be found in the appendix. 
For example, if we consider two different positive eigen- 
values uox and uoxf , the corresponding eigenvectors satisfy 
the following orthogonality relation 

2N 

{^x^y + c^o')«^ = 0, A / y, (35) 



which is derived in the appendix. The normalization of 
a\ is fixed by the following commutation relation 



(36) 



-UJX' 



^UJx 

hx' 



which has been simplified by using Eqs. (34) and (35) 



We are interested in representing the operator (j^ in 
terms of the raising and lowering operators ax and so 
we can relate the displacement dR^ in the plane to the 



collective mode operators ax and a\. Using Eq. (|30|), 



is related to ax and a\ via 



= 2i 



mui, 



1/2 



h 



(37) 
(38) 



in which the relation = imujQ^a\/uox^ Eq. (A6) and 
Eq. (A7) are used. Hence, 



X 



(39) 



This expression is crucial for the discussion of the planar 
modes in a quantum simulation. The (3 component of the 
ion displacement associated with phonon coherent state 
in Heisenberg picture can be evaluated as 



A/3v 



2% 



E 

A:WA>0 



(40) 



= -2n |</'A||«f |sin(wAt + <5A) (41) 

A:u;a>0 

where is the complex eigenvalue for the coherent state 
\(j)) (a A 10) = 0a|0)) 7 is the phase angle for mode A, 
and \(j)x\ is the average phonon occupation for mode A, 



a 



/3A 



Xli/ ^x^j^ i^ projection of A state at site j 



along the (3 orientation. 



D. Effective spin Hamiltonian with axial phonon 
modes 



We discuss the generation of effective spin-spin cou- 
plings from an optical dipole force within the context 
of the Penning trap experiment at NIST, which uses 
ions of Beryllium (modifications for other systems would 
be straightforward, but are omitted here). Consider 
the hyperfine qubit states |tz) = | ^81/2, = 1/2) and 
|lz) = | mj = —1/2) in a ^Be^ ion. The qubit level 

splitting in the presence of a magnetic field of magnitude 
Bz = 4.46 Tesla (due to Zeeman effects) is approximately 
27r X 124 GHz. Even though the ^Be^ nuclear spin is not 
a spin singlet, the nuclear spin dynamics can be com- 
pletely frozen out by optically pumping the nuclear spin 
to the single nuclear spin state mj = -1-3/2 throughout 
the duration of an experiment [H [181 HI] , which is what 
we assume is done here. 

A spin-dependent force along the Z axis of the Bloch 

which is im- 



sphere is generated through a a gate 
plemented by a coupling of the qubit states with the ex- 
cited states (Itz) ^1 ^P3/2, rnj) and \iz) ^\ ^P3/2, rnj)) 
via two off-resonant laser beams with different angular 
frequencies uju •, and wavevectors ku , Icl respectively. 
The AC Stark shift due to the interference of two linear- 
polarized beams for an ion labeled by j generates the 
following spin-dependent force (after adjusting the po- 
larization of the laser beams) [8 : 



Hod 



Fo 



N 

^ 5k, 



cos{5kzj — lilt) a f 



(42) 



where the transverse wavenumber |ku — IclI = Sk, ~ 
2/c sin(^i^/2)(|/c[/| = k ~ |^l|) is determined by the 
magnitude of the wavevector difference between the two 
beams, 6r/2 is the angle between the beam with respect 
to the tangential orientation of the two-dimensional crys- 
tal plane, /i = cjf/ — cj^ is the beatnote frequency given by 
the frequency difference of the two beams, and zj = SRj 
is the displacement of the transverse phonons along the 
axial direction z. The Pauli spin operator shows that the 
force is equal in magnitude, but opposite in direction for 
each of the two spin states of the Beryllium ion (hence it 
is known as a spin-dependent force). In the Lamb-Dicke 
limit, we have \SkzSRj \ <C 1, so the interaction Hqd can 
be further reduced to the following form 



N 



Hod = - ^ FoSR] sm{fit)af 



(43) 



after dropping a term with no SRj dependence. In this 
limit, the spin-dependent force —Fq sin fitaj is spatially 
uniform. 

The effective spin Hamiltonian and the time depen- 
dent evolution of the wave function from the Hamiltonian 
Hod have been rigorously studied when the system is 
cooled to low temperature to start from the phonon vac- 
uum ^21 [30] immediately before the onset of the quantum 
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simulation. The evolution of the entangled spin-phonon 
states are captured by the time evolution operator 

U{t,to) = eM-iHlh{t-to)/n] eiiv[-iWi{t)in]U,pin{t,to) 

(44) 

with the operator Wi{t) determined by Wiit) = 
j^^dt'Vi{t') where the interaction Vi{t') is the optical 
dipole interaction term expressed in the Heisenberg pic- 
ture of the bare phonon states (also called the interaction 
picture) via 

Vi{t) = exp[in^^{t-to)/h]noDF{t)exp[-in^^{t-t^^ 

(45) 

in which 7i^^ is the (time-independent) bare planar 



phonon Hamiltonian given in Eq. (20). As a consequence, 
we find the effective spin Hamiltonian Uspin{t^to) is dic- 
tated by an Ising spin Hamiltonian 



Uspinit^to) = Tteyip 



N 



to 



dt' Jjr(t>f4 



(46) 



with time-dependent exchange integrals. 

These spin-spin exchange integrals are given by [13] 



F2 ^ l)Zui^zu 



J J' 



4m ^ — ' a 



X [l + cos2/it- 



2/i . 



sincj^^ytsin/it]. (47) 



the errors for a quantum simulation mediated through 
the axial phonon modes when the momentum difference 
fi(ku — Icl) of the two Raman beams is not oriented pre- 
cisely along the z axis due to laser alignment errors. Nev- 
ertheless, one can also involve the planar modes on pur- 
pose for quantum simulation, which may be useful when 
the planar modes are far away from the axial modes in 
energy. To simplify our discussion, let us choose the ori- 
entation of the planar spin-dependent forces due to the 
off-resonant laser beams to be along the x-axis in the 
laboratory frame. 

Analogous to the earlier discussion for axial phonon 
modes, the AC Stark shift along the x-axis for the planar 
phonon modes is associated with the effective momentum 
transfer fi5kxX due to the photons via 



N 



OD 



^t)af. (48) 



Notice that the ion coordinates xj {t) = x^ {t) + 6Rj {t) in 
the laboratory frame are determined by the rigid-body 
rotation of the ion crystals with the x axis chosen along 
the orientation of the effective wavevector ku — kL = 
SkxX in the laboratory frame and fluctuations 5r^(t) in 
the laboratory frame. Hence, in the Lamb-Dicke limit 
\5kx5r^\ <C 1, the optical dipole interaction Hqd is de- 
scribed by 



N 



Hod = Y,FoSr^t) sm{Sk,x''^{t) 



(49) 



The detailed phonon mode properties cj^zy, ^j^, and b^y 
determine the values of the different spin-spin interac- 
tions. The sign of the effective spin-spin interaction 
Jjj'{t) depends on the redness or blueness of the laser 
detuning ji with respect to each phonon mode and the 
range of the interaction can be tuned depending upon the 
magnitude of the laser detuning away from all the axial 
phonon modes. In our numerical discussion, we will show 
these effects with direct numerical calculations. 

Note that the analysis becomes more complicated if 
there is a transverse magnetic field present in addition to 
the spin-dependent optical dipole force. This is needed 
for quantum simulations that employ adiabatic state cre- 
ation. Nevertheless, we do not discuss the evolution of 
the system further here when there is a transverse mag- 
netic field present. Details of the case of the Paul trap 
can be found in Ref. [3Ql 



where the phase SkxX^{t) — fit for the planar modes is 
modulated by the rotation of the crystal for the ion 
coordinates in the laboratory frame given by x^{t) = 
cos{—ujt-\-(j)^) with the static phase (j)^ determined by 
the equilibrium configuration in the rotating frame with 
respect to the orientation of the spin-dependent force act- 
ing along the direction x at the initial time to = 0. Hence, 
the optical dipole interaction in Eq. (49 ) can be expanded 
in harmonics of Q. We expect the effective spin- spin 
interaction due to the planar modes to be sensitive to 
the rotational angular frequencies and the laser detuning 
/i. The value of the spin- spin interaction is sensitive to 
the structure of the different planar modes. Using the 
standard partial wave expansion [31 , we can expand the 
function sm{SkxX^{t) — jat) in harmonics of Q as 

sm{Sk^x°{t) - fit) = J2fiiWkxR°) 



E. Spin-dependent interaction by planar modes 

Let us now discuss the situation where the momen- 
tum transfer from the two laser beams lies in the crystal 
plane. One complication of this set-up is due to the rota- 
tion of the crystal with respect to the laser beams in the 
laboratory frame of reference. Examining the planar- 
mode coupling is also important in order to estimate 



X Pi[cos{-Qt- 



(50) 



where the function fi satisfies fi{t) 
for an even integer and fi{t) = 



i^{2l + 1) sin /it 
-^(2/ + 1) cos /it 
for an odd integer The function ji is the spheri- 
cal Bessel function of the first kind, and the function 
Pi[cos{—Qt-\-(pj)] is the Legendre polynomial in harmon- 
ics of Q. Notice also that the fluctuations 5rj are related 
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to the fluctuations due to phonons in the rotating frame 
as 6rj = dRj cos + dlij sin Qt. It is tempting to sim- 
phfy the discussion further based on the behef that one 
can tune the laser such that only one targeted planar 
mode is closest in energy. Considering the narrow band- 
width of the planar modes, which is roughly an order 
of magnitude narrower than the bandwidth of the axial 
modes, the optical dipole interaction due to the planar 
phonon modes cannot usually be described by a single 
phonon mode even when the laser detuning is close to 
that mode. Therefore, we do not take this strategy for 
the planar modes. In the following discussion, our deriva- 
tion of the effective spin Hamiltonian is not restricted to 
any particular laser detuning or phonon band. 

Following the procedure we mentioned earlier, the ef- 
fective spin Hamiltonian due to the optical dipole inter- 
action can be extracted by the commutation relation as 



dipole interaction is given by 

N 

Hod = - ^FqSR^ sm{Sk^R^j cos - fit)af , (54) 



in which the relation x^{t) ~ coscj)^ is used. By direct 
evaluation of the effective spin Hamiltonian in Eq. (51), 



Ha 



2h 



(51) 



which 



the 



deflnition 



Vj{t) 



X fi{t)Pi[cos{-nt ■ 



(52) 



the effective Ising Hamiltonian Hspin can be cast into the 
form 

Hspin = riFl Im[ara^7/(i)sin(4-Aii)e-*'^^*" 

jj'Xvv' 



VXVLXV ^Zj ^Zj 



//(t) = f dt'e'"''' sin{^° - fit'). 
Jo 



(55) 



exp[iH^^t/h]HoD{t) exp[-iH^^t/h] and Wi{t) = 

j^dt'Vi{t') are applied again. After tedious algebra, 
we arrive at the following expression for the effective 
Ising spin Hamiltonian Hspin = ^jj' with 
spin-spin interactions given by 



Since the contribution from fast oscillatory frequency 
components in Hspin is expected to be small because the 
contribution to the quantum phases are averaged out, we 
disregard the high frequency modes ±2/i under the rotat- 
ing wave approximation. We arrive at the main conclu- 
sion that the effective spin simulator in the slow limit 
{QT <C 1) can be described by the static Ising model 

with the interaction Hspin = ^jj' '^jj'^j^f which 
the static Ising coupling is given by the expression 



7° 



E 



ujx COS (pjjfKela^^a^ } 



Xvv' 

lism.(j)jj>l-m{a*^ a\ } 



WW 



(56) 



in which the indices /3, {3' run through y components 
of the normal modes and the time dependent functions 
h^i , g^^i are deflned by the following expressions 

g^'l = cos(m) / dt7^(tOcos(mOe'^^*'Pz[cos(-m' + (^^^)] 
Jo 

= sin(m) / dt'fi{t')cos{nt')e'^^^'Pi[cos{-nt'^(l)^j^ 
Jo 

h^^i = cos{nt) [ dt7^(tOsin(mOe^^^^'Pz[cos(-^]t' + (^^°)] 
Jo 

= sin(m) / dt7z(tOsin(mOe'^^''Pjcos(-m' + ^^°)]. 
Jo 

(53) 

It is possible that results drawn from these general for- 
mulas could be examined in benchmarking experiments 
such as the measurement of the spin-echo response or an 
average spin-spin interaction, when the laser beams have 
an effective wavevector parallel to the crystal plane. 

Another interesting limit is the occasion in which the 
simulation time T is much shorter than the inverse of the 
rotational frequency QT <C 1. In this limit, the optical 



and the angle = (j)^ — (j)^, is determined by the ion 
equilibrium configuration. Because the function sinOjj' 
is odd for the permutation j ^ j' and the expression 
^vv' WW ^^^^^ permutation v ^ v' \s unchanged, 
only the first term in J^-, contributes to the effective Ising 

Hamiltonian Hspin after the summation over the indices j 
and j' . We therefore arrive at the main conclusion for the 
planar mode spin Hamiltonian Hspin = ^jj' 'Jjj'^f^j' 
slow rotation wT <C 1 where the spin-spin interaction is 
described by 



7^ 



(57) 



The effects of the distribution of the ion equilibrium po- 
sitions from the interaction in Eq. (54) can be inter- 



preted as a varying phase fiuctuation due to the exter- 
nal symmetry breaking from the laser beams. In our 
later numerical discussion, we show two planar mode 
branches. One branch has much lower energy than the 
other branch. The quantum simulation with phonons 
in the lower branch is particularly relevant for exper- 
imental realization when the laser is closely detuned 
from that branch. In this case, the larger spin-spin ex- 
change interaction (which is inversely proportional to 



10 



11^ — uj\^ 2ujx X — ^a) where ji ^ uox guarantees that 
the simulation can be reahzed at a much shorter time 
scale before the decoherence from the environment takes 
place. 



III. NUMERICAL RESULTS 

We have discussed above the theoretical formulation 
for the description of phonon-mediated quantum spin dy- 
namics. Due to the inhomogeneity and finiteness of the 
experimental systems, numerical studies are necessary to 
understand the behavior and facilitate a detailed under- 
standing of current and future cold-ion experiments in 
Penning traps. 



A. Typical system parameters 

We consider ^Be^ ions trapped in a static magnetic 
field with a rotating wall potential yvt^r^cos2(^ + Vtt) 
{VL > 0). The axial trapping frequency due to the end- 
cap potentials eVo = is fixed at the value uj^ = 
27r X 795 kHz (in units of rad/s) as typically applied in 
experiments. We measure all angular frequencies in units 
of ujz- For example, we define the frequency scale ujw 
associated with the rotating wall potential Vw by the 
relation ujw = \/2e|lV|/m. We choose cases with very 
weak, weak, and strong rotating wall potentials given 
by the corresponding values uw = O.Ola;^, 0.04^;;^, and 
0.07uJz- The cyclotron frequency ujc associated with the 
magnetic field is uJc = eBz/m. With the magnetic field 
Bz = 4.5 Tesla, the value of the cyclotron frequency is 
given by uJc = 9.645co'2. The Beryllium ion has an atomic 
mass m = 9.012182u, where u is the atomic mass unit, 
and has a positive unit charge e = 1.60217646 x 10~^^ 
Coulomb. 

Based on our theoretical discussion above, the system 
is energetically stable when the rotational frequency Q 
due to the rotating wall potential lies between the de- 
confinement transition frequency uj^c and the one-to-two 
plane transition frequency uji2- The deconfinement fre- 
quency ujdc is determined by the criterion /^a = 0. This 
deconfinement criterion can be recast into the following 
form: 

^dc = Y-\lf-f-<- (58) 

The one-to-two plane transition frequency uji2 needs to 
be reliably determined numerically by the instability of 
the axial phonon frequencies. We can map out the one- 
to-two plane instability by calculating when the axial 
phonon frequencies first become imaginary, as shown in 
Fig. [2] for Vw = 0. The one-to-two plane transition fre- 
quency depends only weakly on the wall potential via the 
change of the equilibrium coordinates due to the presence 
of the wall. 




z 

FIG. 2: One-to-two plane transition frequency for different 
numbers of ions N. The vertical dashed line shows the de- 
confinement frequency for Vw = 0. For larger values of W, 
the vertical line moves to the right, reducing the frequency 
range for the rotating wall, where the single plane is stable, 
while the one-to-two plane transition line hardly changes. 



In our numerical discussion, we choose three rota- 
tion frequencies Q = UJ12 — A, = 0.5(6^^^ + ^^12), 
Q = uJdc + ^[^dc depends on Vw] with the frequency 
offset A = 27r X 200 Hz where the frequency 00^^ is the 
deconfinement frequency at zero rotating wall potential 
{Vw = ojw = 0). Equivalently, instead of the rotational 
frequency, we report the effective trapping frequency in 

the crystal plane cjg/ / = \J — — ^ where the re- 
lation eVo = 0.5mcj^ is used. Thus, we can report the cor- 
responding scaled effective frequency uj^ = 0.21ujz for the 
high rotational frequency case 60^12 — A, the corresponding 
mean effective trapping frequency Um = 0.166^^ for the 
mean rotational frequency case 0.5(0;^^ + ^12)^ and the 
low effective trapping frequencies for the low frequency 
case Q = uJdc + A with very weak, weak and strong rotat- 
ing wall potentials as oo^^ = 0.050;^, = O.OGuJz^ and 
cjf = 0.09uJz^ respectively. We have three cases for the 
low frequency case because the deconfinement frequency 
depends strongly on the wall potential near the low fre- 
quency, so it must change with the wall potential; for the 
other two cases the dependence is so small that we ignore 
it. 



B. Equilibrium configurations 

To seek for the equilibrium configuration of crystal- 
lized ions in the trap, we generate a seed lattice of finite 
size, which takes the form of a triangular lattice near 
the center and has the boundary as a closed hexagonal 
shell, with six facets. We look for a stable equilibrium 
configuration which minimizes the Hamiltonian in the ro- 
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tating frame, assuming none of the ions are moving. We 
find this procedure always guarantees that we find a lo- 
cal minimum that looks like a triangular lattice in the 
center, just as is found in experiment but has an ellipti- 
cal boundary which is shaped by the equipotential of the 
effective trapping potential (this approach appears to be 
similar to that in Ref. 16 , but we use a different potential 
by including the rotating wall, while it is different from 
the approach used in Ref. [ITj). The minimization pro- 
cedure employs both a locally calculated gradient and a 
locally calculated Hessian matrix, similar to the method 
described in Ref. [32] We used the MatLab optimiza- 
tion toolbox and our own code to calculate these minimal 
configurations, and results between the two codes agreed 
very well. 

The construction of the seed lattice is based on a closed 
shell construction analogous in spirit to closed-shell elec- 
tron configurations for atoms. More specifically, for any 
particle number such as N = 217, we first determine how 
many closed hexagonal shells can be created with that 
number of ions. The lattice is a closed hexagonal lattice 
if the boundary of the initial seed lattice forms a perfect 
hexagonal shape. The number of the closed hexagonal 
shells S is related to the number of ions N for the seed 
lattice by 



S = 



V9 + 12(7V-1) 



(59) 
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3: Spatial configuration for a seed lattice with = 217 
The spatial dimension is scaled by the characteristic 



Vo 



V ™z J 



length in the axial direction z as I 

The red cross denotes the point where the external potential 
due to the electrodes vanishes. The minimal ion configura- 
tion often, but not always, has an ion located at this energy 
minimum. 



If the number of ions N cannot fill an integer number 
hexagonal shells, the outermost ions are placed on an i 
complete outer hexagonal ring according to the minim 
potential energy at each of the outer ring sites (the ene 
gies are different because the outer ring is not elliptic 
in shape, especially when Vw is nonzero). For the ca 
with = 217, an integer number of closed hexagon 
shells are generated in the seed lattice {S = 8) which 
shown in Fig. [3j 

In Fig. |4j we show spatial equilibrium configurations 
different rotational frequencies with the same fixed wei 
rotating potential uw = 0.04cj2. We expect that t] 
anisotropy introduced by the rotating wall potential u 
is minimal when the effective trapping potential domi- 
nates over the rotating wall potential oow As shown in 
the Fig. [4jc), we clearly observe that the ions are ar- 
ranged much more isotropically with a much higher den- 
sity than Fig.|4ja) and Fig.[4]^b) due to the large effective 
trapping in the plane. Note how the seed lattice smoothly 
evolves from a nice triangular lattice in the center to a 
boundary that takes the minimal, elliptical shape, of the 
potential at the edges, in the stable equilibrium positions. 

To further quantify the change in isotropy, we define 
the distortion of the crystal as the aspect ratio of the 
major axis to the minor axis of the elliptical edge. Fig. |5] 
shows the dependence of this distortion on the rotational 
frequencies oOeff and rotating wall potential ujw With 
a fixed rotating wall potential, we observe the lattice is 
more distorted with a weaker effective frequency oJeff- 
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FIG. 4: Equilibrium structures for different rotating frequen- 
cies: (a) Low effective trapping frequency u]^ — O.OGcj^. (b) 



Medium effective trapping frequency ujm 
effective trapping frequency ujh = 0.21ujz 
is uw = 0. 04:00 z for all cases. 



= 0.16cj^. (c) High 
The wall potential 



The effect is much larger for the strong rotating wall 
potential uw This shows that the role of the rotating 
wall potential is not always negligible. Its effect can only 
be neglected when the ratio between the potential and 
effective trapping frequency uw/^ef f is smaller than one. 

Fig. [6] shows how the ion spacing of the equilibrium 
structure changes as a function of the distance from the 
origin. To compute the lattice spacing for a given ion, we 
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FIG. 5: (Color online.) Dependence of the crystal distor- 
tion ratio (or aspect ratio) on the rotational frequency ujeff 
and the rotating wall potential cuw- The black-solid curve 
represents a strong rotating wall with uw = O.OTcj^. The 
green-dashed curve is a weak rotating wall potential with 
uow = 0.04:UJz. The red-dotted curve is a very weak rotat- 
ing wall potential with ujw = O.Olcj^. The colored dots are 
the principle cases we examine with more detailed calculations 
for much of the remaining numerical work. 
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FIG. 6: Ion spacing dependence as a function of the dis- 
tance R from the origin. The black, green, and the red 
symbols denote respectively the data with a strong rotat- 
ing wall potential uJw = O.OTcj^, a weak rotating wall po- 
tential uw = 0.04:UJz, and a very weak rotating wall potential 
(jJw = O.Olcj^. The values of the corresponding effective trap- 
ping frequencies cje// are shown in the legend. The same 
parameter sets as in Fig. |5] are used. 



find its nearest neighbors using a Delaunay triangulation 
algorithm and then calculate the mean distance from that 
ion to those nearest neighbors. Each data point repre- 
sents the mean nearest-neighbor ion spacing for each ion 
at a different radius R. When we have a strong rotating 
wall and a weak effective potential, we expect the ions 
are more widely spread out and are pushed more along 
the weakly confined direction. This can be seen in the 
data sets with colored triangles, which have a larger mean 
nearest-neighbor distance and a larger spread about the 
average. For the data sets with colored squares and cir- 
cles (weaker rotating wall potential), we observe much 
smaller mean nearest-neighbor distances with much less 
spread. 



C. Normal modes 

In this subsection, we discuss the nature of the nor- 
mal modes including the energy spectrum and properties 
of the eigenvectors. Here we present the eigenfrequen- 
cies from each of the three phonon branches (one branch 
from the axial degrees of freedom and two branches from 
the planar degrees of freedom). We find in general that 
the three branches are well separated in energies with the 
bandwidth strongly dependent on the equilibrium config- 
uration of the ions. The upper planar modes (shown in 
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FIG. 7: (Color online.) Eigenfrequencies of the upper 
(cyclotron-like) branch of the planar modes. The same pa- 
rameter sets and notation are used as in Fig. [5] 



Fig. [7|) have the character of the high frequency cyclotron 
motion renormalized by the Coulomb interaction between 
ions and the lower planar modes (shown in Fig. [9| have 
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the character of the low frequency magnetron motion. 
The eigenfrequencies of the axial mode branch shown in 
(Fig. [8| lie in between the two planar mode branches. In 
each figure, the modes are sorted by ascending frequency. 
In general, eigenfrequencies depend mostly on the effec- 
tive frequency cjg// and are less sensitive to the small 
perturbation due to the rotating wall potential uow even 
though the character of the eigenvectors can be quite 
sensitive to the rotating wall potential. 
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FIG. 8: (Color online.) Eigenfrequencies of the axial modes. 
The same parameter sets and notation are used as in Fig. [5] 

In Fig. [7| the dependence of the cyclotron-like eigen- 
frequencies on the rotating wall potential and rotational 
frequency is shown for the nine typical cases we will be ex- 
amining in our numerical results. Cases are represented 
by the same color for the same ujw with the same sym- 
bol shape for the same cjg//. Note how it is the effec- 
tive frequency that primarily determines the frequency 
spectrum, with a separation due to the rotating wall po- 
tential only becoming clearly visible when that potential 
is strong. In addition, note the lowering of the lowest 
eigenfrequency with large cjg//, which signals a collective 
feature due to Coulomb interactions (colored circles). 

In Fig. [8) we observe one axial mode (the highest one) 
which has a universal eigenfrequency with uj = uOz inde- 
pendent of the parameter sets. This is due to the fact 
that the uniform center of mass motion along the axial 
direction does not cost any Coulomb energy and has the 
same trapping energy in the axial direction so it is inde- 
pendent of the rotating wall potential applied in the crys- 
tal plane. This is shown as the intersection of the highest 
axial modes for all the parameter sets at the highest axial 
frequency {uj = uOz)- The other phonon modes have lower 
eigenfrequencies than the center of mass (COM) mode. 
This can be interpreted as the reduction of the Coulomb 



interaction due to the increase of the average ion spac- 
ing from the phonon displacement when the wavevector is 
nonzero (for the COM mode, the phonon does not change 
the relative displacement of the ions, so the Coulomb in- 
teraction is unchanged). This interpretation also agrees 
with the observation that the axial mode branches for 
large oo^ff lie lower in energy as most explicitly shown 
with the colored circle data. The phonon frequencies near 
the COM mode have been verified against experimental 
measurements using a technique that can measure the 
effective temperature per phonon mode and the actual 
frequency of each phonon mode [21 . 

In Fig. [oj we plot the lowest (magnetron-like) planar 
modes for the same nine cases. When there is no rotat- 
ing wall potential, the lowest-frequency mode is that of 
a rigid rotation, which costs no energy since the poten- 
tial in the plane depends only on the radius in that case. 
This zero energy mode is the Goldstone mode due to the 
spontaneous rotational symmetry breaking of the equilib- 
rium configuration. As we break the rotational symmetry 
by turning on a wall potential, the Goldstone mode "ac- 
quires mass" and becomes nonzero, but often remains a 
very low frequency mode with a near rigid rotation char- 
acter. The bandwidth of this branch of planar modes is 
quite sensitive to the effective trapping frequency ujeff 
as can be shown by examining the the colored circle and 
colored square cases. The bandwidth is only slightly de- 
pendent on the rotating wall potential as can be seen 
with the colored triangular cases. 

As the rotational frequency is increased relative to 
the axial frequency, the confinement in the plane be- 
comes tighter and tighter until, suddenly, the system 
becomes stabilized in a two-plane configuration. As we 
approach this critical frequency, the axial mode phonon 
bandwidth grows tremendously, as the lowest axial fre- 
quency marches down and eventually becomes lower than 
the maximal magnetron-like planar phonon frequency 
(which increases with the rotational frequency but not 
as rapidly), as shown in Fig. 10 The overlap occurs (at 



uJeff = 0.21896^^) quite close to the critical one-to- two- 
plane transition {oJeff = 0.22020;^) for a rotating wall 
potential corresponding to ujw = O.OAujz- 

We now discuss the nature of the eigenvectors for the 
different phonon modes. In general, we find phonon 
modes to either be collective, where all of the ions move 
with a long- wavelength, similar to the drumhead modes 
of the continuum, or to be localized, where the motion 
is confined primarily to a small fraction of the ions. It is 
the long-wavelength modes that are most important for 
quantum simulation, so we focus on them. For the axial 
branch, these modes lie close to the COM mode (at high 
frequency), while for the magnetron-like branch, they are 
close to the rigid rotation mode, at low frequency. In 
Fig. [TT] we show the four highest frequency modes of the 
axial branch and in Fig. [T2j we show the two lowest fre- 
quency modes of the upper (cyclotron-like) branch and 
lower (magnetron-like) branch planar modes, all for the 
same parameters as in Fig. [4j Unlike Fig. [4j our focus 
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is on the character of the eigenvectors instead of the ab- 
solute size of the equihbrium configurations, so we have 
rescaled the positions to make the ion positions have sim- 
ilar widths in all panels. Similarly, the color scales used 
for the axial modes and the arrow lengths used for the 
planar modes are both adjusted for each panel to high- 
light the details of the motion and should not be used 
to compare relative ion displacements from one panel to 
another. 

The leftmost 



We start with the axial modes in Fig. 11 



mode is the COM mode which involves uniform motion 
of all of the ions for all cases. If we start examining 
the next two modes in Figs. pT]^j)-(k) and Figs, [lljf)- 
(g), we see that they are two nearly degenerate states 
due to the weak anisotropy generated by the rotating 
wall, and they are almost related by a rotation by 7r/2. 
Interestingly, as the effective frequency gets smaller in 
Figs. 11 ^b) and (d), we find that the anisotropy forces 



another mode in Fig. [TT[c) in between these two closely 
related modes. The remaining mode in Figs.pTj^h) and (1) 
is close to the next drumhead mode which has two lines of 
nodes. This feature is lost with strong enough anisotropy 
in Fig. 11 'c). The anisotropy increases because the slower 



rotating crystals feel the effect of the wall potential more 
strongly. These first few modes agree well with a contin- 
uum elastic theory for the normal modes . One 
interesting feature of all the eigenstates is that the spatial 
variations are periodic along the circumference of the ion 
cloud but there is no constraint along the radial direction. 
The quantum dynamics of any axial normal mode v can 
be understood as the collective quantized oscillation with 
the eigenfrequency uj^y along the z axis with the given 
snapshot of the eigenvector shown in Fig. 
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FIG. 9: (Color online.) Eigenfrequencies of the lower 
(magnetron- like) branch of the planar modes. The same pa- 
rameter sets and notation are used as in Fig. [5] 




FIG. 10: (Color online.) Comparison of the phonon band- 
widths for the axial (red) and lower (magnetron-like) planar 
(blue) branches for a fixed weak rotating wall potential. As 
the effective trap frequency increases, the lowest frequency of 
the axial mode decreases as the largest (magnetron-like) lower 
branch planar frequency increases until they overlap very near 
the one-to-two plane transition. The uppermost axial mode 
is fixed, while the lowest planar mode does not change signif- 
icantly with CJe//- 

Due to the presence of the magnetic field, the ion dis- 
placement in the xy plane is chiral arising from the cou- 
pling of the dynamics between the x and y components of 
the same ion. We show that the lower branch of the pla- 
nar modes is the slow collective renormalized magnetron 
modes inherited from the single ion features under E x B 
drift. The higher branch of the planar modes; however, 
are the fast renormalized cyclotron modes mostly inher- 
ited from the the single ion features under the the strong 
static magnetic field B^z. The snapshot of the lowest 
two eigenmodes in the two branches are shown in Fig. [12) 
The corresponding magnitude of the displacement for a 
mode A is represented by the length of the arrow and the 
direction of the arrow represents the displacement of the 
ion at the time when the snapshot is taken. 

However, to understand the character for the pla- 
nar mode A, it is revealing to show the time evolu- 
tion of the snapshots based on the displacement of the 
phonon coherent state due to the mode A projected 
on the ion lattice positions {{6 (t)) , {S Rf (t))) ex 
{\a^ ^\ cos ujxt,\a^^\ cos ouxt) where the phase Sx can be 



chosen to be — 7r/2 in Eq. (41) for our purpose. We pro- 
vide a movie of this in the supplementary material [35]. 
Here we summarize what we have observed for the dy- 
namics of the planar modes. Figs. 12 ^b), (f) and (j) 
show the lowest collective cyclotron modes for different 
effective trapping frequencies. We observe that each ion 
shows clockwise cyclotron motion around its equilibrium 
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FIG. 11: (Color online.) Four highest frequency axial eigen- 
vectors for three different effective trapping frequencies with 
a weak rotating wall potential uw = 0.04cc;^. Subplots (a)- 
(d) present the case with a low effective trapping frequency 
uj^ = 0.06(jJz, (e)-(h) present the case with a mean effective 
trapping frequency cjrn = O.lGcj^, and (i)-(l) represent the case 
with a high effective trapping frequency ujh = 0.21(jJz- The 
color in the colorbar represents the scaled value of the nor- 
malized eigenvector b^^ with respect to its maximal positive 
component for the mode 
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FIG. 12: A snapshot of the spatial eigenvectors for the lowest 
two modes in each planar branch at the fixed rotating wall 
potential ujw — 0.04c(;2. Subplots (a)-(b), (e)-(f), and (i)-(j) 
show the lowest two planar modes in the upper branch with 
the effective trapping frequencies uoY — O.OScj^, cjm = O.OGcj^^, 
and uJh — O.OQcj^ respectively. Subplots (c)-(d), (g)-(h), and 
(k)-(l) on the right panel are the corresponding data for the 
lower branch of the planar phonon modes. 



position in the rotating frame. The second lowest cy- 
clotron modes shown in Figs. [12] (a), (e), and (i) also show 
clockwise cyclotron motions for each ion around its equi- 
librium but with nonuniform fluctuations on the mag- 
nitude of the displacement. The quantum dynamics for 
the renormalized cyclotron modes are described by the 
dynamics of the quantized cyclotron orbits. Since these 
orbits have much higher energy than the lower branch 
of the planar modes, the quantized cyclotron orbits have 
much smaller average radius. 

^d), (h), and (1), we show the global shear 



In Figs. 
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mode (lowest mode in the lower branch), which is adi- 
abatically connected to the global rotational mode at 
very weak rotating wall potential, with zero frequency 
in the rotating frame. The normal mode displacement in 
principle oscillates back and forth infinitesimally slowly 
along the elliptical shell of the equilibrium configurations. 
The snapshots of those modes are taken at different ini- 
tial conditions and therefore with different orientations of 
the displacement. In Figs.[T2|c), (g), and (k), we observe 
multiple fragmented shear domains with different signs of 
circulation. The circulation with one sign of circulation 
is always suppressed at low energy. 



D. Ising spin-spin interactions 

Now that we have described the phonon properties, we 
are ready to describe the properties of spin-spin interac- 
tions when the laser is detuned close to different branches 
of the phonon modes. 

Let us focus on the axial phonon modes first which have 
been used to benchmark the mean spin-spin interactions 
between hundreds of spins in a Penning trap [8]. We 
calculate the static spin-spin interaction J^j between ions 
i and j based on Eq. (47). 

In Fig. [13) the detunings (5 = fi — Uz range from 
O.OOOlujz to ujz^ all detuned to the blue of the COM mode. 
In this case, all of the interactions are antiferromagnetic 
and they decay with a power law depending on the dis- 



tance between two spins {Jij oc 



' ij J- 



The power law 



increases from a = 0, when the detuning is very close 
to the COM mode to a = 3, when the detuning is far 
away from the COM mode. We observe that the case of 
weak effective trapping allows for a much more rapid ap- 
proach to the limit of dipole-dipole interactions {a = 3) 
at large detuning, as shown by the data sets with colored 
triangles. 

We next examine what happens when we detune the 
spin-dependent optical dipole force to lie in between dif- 
ferent axial modes. In this case, the interaction Jij 
does not reveal well-defined power law behaviors and is 
frustrated with different signs; therefore, we report his- 
tograms of the spin-spin interactions. We choose the fol- 
lowing four detunings: the mean of the first and second 
lowest frequency axial modes, the mean of the 72th and 
73th axial modes, the mean of the the 144th and 145th 
axial modes, and the mean of the first and second highest 
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1500 





FIG. 13: (Color online.) Detuning above the axial modes 
(blue of the COM mode). As in Fig. 2(b) of Ref.EI a detuning 
= 6-\-uJz above ujz yields a power- law- like decay of Jij oc r~-^ 
as a function of the ion distance, Vij. The legend presents 
the values for the effective trapping frequencies in units of 
ujz. Here we plot the exponent of the fitted power law as a 
function of the detuning. 



frequency axial modes. In Fig. [T4| we observe the distri- 
bution of the interaction J^j is more symmetric about 
zero when the detuning /i is located near the central part 
of the axial modes such as cases (b) and (c). In addition, 
by our more detailed studies, we found that the spin-spin 
interaction Jij is randomly disordered with distance but 
is correlated with the polar angle between ions. This 
can be interesting for further studies of disordered spin 
dynamics with the spin-spin interaction. However, the 
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FIG. 14: (Color online.) Detuning inside the axial mode 
branch with an effective trapping ujrn = 0.1 Gcj^ and the weak 
rotating wall potential uw = — O.OGcj^ for N — 217. (a) 
Detuning 6 = 0.00567ujz below uJz exactly between the two 
highest modes with 99.93% of the Jij in the plot, (b) Detun- 
ing S — ^. 117 Aujz below UJz exactly between the [2A/'/3j and 
[2iV/3j + 1 mode with 99.28% of the Jij in the plot, (c) De- 
tuning 5 — ^.1921ujz below uJz exactly between the [A/'/3j and 
[Ar/3J + 1 mode with 59.73% of the Jij in the plot, (d) De- 
tuning 5 — 0.292Qujz below uJz exactly between the two lowest 
modes with 63.65% of the Jij in the plot. The inserted green 
curves are generated from Fig. [S] with the arrow indicating 
where the detuning lies. 



nar mode at /i = 1 x lO^^cj^, in between the lower and 
upper branch planar modes at /i = cj^, and above the 
upper branch planar modes at /i = l^uOz^ corresponds to 
Figs, p^a), (b), and (c), respectively. The darker colors 
(red, green, blue) represent positive Jij and the lighter 



interaction Jij becomes asymmetric when the detuning colors (magenta, yellow, cyan) represent negative Jij val- 



II is located elsewhere. We also notice that the fluctua- 
tions of Jij are larger when the detuning is closer to the 
high frequency portion of the axial modes. 

Below the axial modes, the detunings are chosen with 
an equal spacing from just below the lowest frequency 
axial mode, to O-Scj^ which is sufficiently far from the 
lower branch planar modes for the effective trapping uOm- 
For large detuning we still observe the majority of 
the spin-spin interaction Jij follows the power law scaling 

Jij (X 1/' 



ues for the same parameter set. One can observe all Jij 
are almost equally distributed among positive and neg- 
ative values for all the cases without showing any clear 
signature of power law behavior. 

In Fig. [T7| we examine whether there is any peculiar 
spatial dependence to the Jij patterns for the case with 
/i = 1 X 1{)~^ijOz' In Figs. 17 'a) and (c), we select the cen- 



r'^j except for the case with very small detuning 



such as the case with 8 = —0.29277uJz in Fig. 15 



Similar to calculating the spin-spin interaction Jij for 
the axial modes, we use Eq. ([57| ) to calculate the Jij 
for the planar modes. In Fig. |16[ we present the planar 
spin-spin interaction Jij as a function of the distance r 
between ions. We show that the couplings Jij do not fol- 
low a well-defined power law scaling, but the distribution 
better resembles the symmetric histogram in Fig. [l4]^b). 
Detuning the lasers near the lowest lower branch pla- 



ter and an outer ion as the reference point respectively 
(shown with the black star ^) and calculate the interac- 
tion Ji^ between each ion i with respect to the reference 
ion. In order to study the dependence of Ji^ on the ion 
distance r, we define the subshells of the crystal as ions 
that have a similar radius from the reference ion. In 
Figs.[l7|;a) and (c), ions are in the same subshell if they 
are illustrated with the same shape and color. From this 
new origin, we define Oi as the polar angle with respect 
to the reference ion and plot the Ji^ from each subshell 
as a function of Oi as shown in Figs, p^b) and (d). Each 
of these figures show that the magnitude \Ji^\ in each 
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FIG. 15: Dependence of the spin-spin interaction as a function 
of the distance r between the ions for a detuning below the 
axial modes (red of the lowest axial mode). The effective 
trapping Um = 0.1 Gcj^ and the weak rotating wall potential 
LUw = (^i^ — 0.06ujz are used. 
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FIG. 16: Dependence of planar spin-spin interaction as a func- 
tion of the distance r between ions for various detunings: (a) 
/i = 1 X 10~^(jJz;(h) ji — uJz'Xc) ji — lOuJz- The darker colors 
illustrate positive Jij and the lighter colors illustrate nega- 
tive Jij. The effective trapping ujm = O.lGcj^ and the weak 



rotating potential ujw 



■ O.OQuz are used. 



shell increases with radius as shown by the largest am- 
plitude of the oscillations for the ions in the same shell 
(data points with the same color data) and Ji^ alternate 



in sign at a certain fixed orientation Oi. In addition, the 
interaction J^^ alternates within the same subshell, which 
is very different from the axial modes. These oscillatory 
correlations appear to exist for other detunings such as 
fi = uoz or (1 = lOcj^ as well. Noticeably, as shown in 
Fig. 17(c), the interaction J^^ is stronger for ions situated 
along the orientation Oi = tt opposite to the orientation 
of the laser wavevector Skx{Oi = 0) due to the symme- 
try breaking of the laser beam. So, in general, the spin 
model realized by coupling to the planar modes creates 
frustration due to randomness and long-range interaction 
for spins, which can be potentially be used as a platform 
to explore spin-glass physics, especially since hundreds of 
ions are accessible in the Penning trap systems. 
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FIG. 17: Angular correlation of Jij with /i = 1 x lO^^a;^. 
(a) and (c) The subshell definition for the ion at the center 
and for an outer ion. (b) and (d) The spin-spin interaction 
Jij of these subshells as a function of radius r. The effective 
trapping uom = O.lGuz and the weak rotating wall potential 
uJw — — O.OGcj^ are used. 



IV. CONCLUSION AND DISCUSSION 

In this work, we have presented a thorough treatment 
of the theoretical background for further development of 
cold ion trap computation in a Penning trap, bringing the 
theory to a similar level to what has been available for lin- 
ear Paul trap simulators for years [131 [M]. particular, 
we have shown how to determine the equilibrium posi- 
tions efficiently, how to find the phonon modes for motion 
about the equilibrium positions, and how to couple the 
phonon motion to spin-dependent forces to determine the 
effective spin exchange between two different ions. Our 
work focused on both axial and planar modes. The latter 
are more involved because they have velocity-dependent 
forces in the rotating frame due to the Coriolis force and 
Lorentz force, and require more complicated methods to 
perform the normal- mode analysis. 

Our results show, that in some cases the system is 
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described by relatively simple Ising models, which can 
have a transverse magnetic field added in order to inves- 
tigate adiabatic state preparation. But, in many cases, 
the spin-exchange patterns are quite complex, and will 
require significant analysis to understand their behavior. 
We were forced to restrict ourselves to examine just a 
handful of different cases in this work, because the pa- 
rameter space is so huge. There may very well be a num- 
ber of interesting new results that can be found by ap- 
plying our approach to different situations than what we 
examined. In any case, this formalism will be quite im- 
portant in analyzing the next generation of experiments 
on Penning-trap-based cold ion computation, especially 
for cases where one maps onto Ising-like models. 

Fully understanding actual experiments will, of course, 
require much deeper analysis than what was given here. 
There are intrinsic errors with respect to the target spin 
Hamiltonian, noise due to external environments, and 
spontaneous emission effects, which may deteriorate the 
quantum simulation. For example, in the presence of a 
transverse magnetic field in an adiabatic state prepara- 
tion protocol, it may not be possible to completely keep 
the phonon degrees of freedom as only virtual excita- 
tions during the simulation [30^. In addition, there are 
errors in the initialization of quantum spin states and of 
the phonon states. For example, when the crystal is not 
sufficiently cold, there will be a background of phonons 
present prior to the start of the simulation. Those can 
modify the fidelity of the quantum simulation and under- 
standing how this works requires future study. 
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Appendix A: Details for solving the quadratic 
eigenvalue problem 

We discuss some of the general properties of the 
quadratic eigenvalue problem here, since this problem is 
not too well known within the physics community. Our 
approach is closely related to the work of Baiko [27 , but 
generalizes his work to the case we need for the Penning 
trap planar phonons. For a more mathematical discus- 
sion, see Ref. 28, 

The planar phonon Hamiltonian (in second quantized 
form) is given by H^^ = Yl'xli ^'^a(<^a^a + ^) with the 
phonon frequencies all positive {ujx > 0). In terms of the 



eigenvectors discussed in Sec. IL C, we can solve the 
collective mode problem by solving the following QEP 
and choosing solutions with positive eigenvalues uJx > 



,^2 { 



0. 



(Al) 



Here A is the label for the different solutions and v is 
the matrix/vector index which is summed over. Hence, 
this is a 2N x 2N quadratic eigenvalue problem, with 
eigenvectors a\. There are 4A^ eigenvalues, and hence AN 
associated eigenvectors. Because the QEP is complex, it 
turns out that the eigenvectors are not orthogonal, and 
hence the behavior is different from a conventional linear 
eigenvalue problem. We will derive the relations that the 
eigenvectors do satisfy by mapping the QEP onto a linear 
eigenvalue problem and using the well-known properties 
of those solutions (this procedure is called linearizing the 
QEP). 

The QEP can be mapped onto a linear hermitian eigen- 
value problem as follows: 



m ^' 







^x(^x 



^\(^x 



(A2) 

The AN X 47V matrix is hermitian since T is real and 
antisymmetric. Note that the QEP, and hence the cor- 
responding linear eigenvalue problem, both have positive 
and negative eigenvalues; for the physical solutions we 
want, we will restrict to positive eigenvalues only. 

It may seem odd that the eigenvalue appears as part 
of the eigenvector for this solution, which might bring 
some confusion about how one solves such a problem. 
If we let Mjyjyf denote the matrix on the left hand side 



of Eq. (A2), and let denote the eigenvector, then the 



equation is a conventional eigenvalue problem Mj^j^f x\ = 
cja^a- After solving this conventional problem, the a\ 
vectors are extracted from the explicit forms given in 
Eq. (A2). One might be concerned about whether the 



eigenvectors can always be written in this form, but this 
always follows from the specific form of the matrix M . 
By taking the complex conjugate of Eq. (Al), we find 



that if a\ is the Ath eigenvector of the QEP with eigen- 
value cja, then a^^ is the eigenvector of the QEP with 
eigenvalue —ujx- Hence the eigenvalues come in ± pairs, 
with corresponding eigenvectors of the linear eigenvalue 
problem {uJxo^\ ^o^a)^ ^^iih. positive eigenvalue cja > 0, 
and {—ujxo^^x ^o^a*)^ ^\ih. negative eigenvalue —uox < 
as long as the eigenvalues uJx and uOq are real. As a con- 
sequence, there are AN eigensolutions (as expected for 
a 4A^ X 4A^ linear eigenvalue problem) and half of them 
have positive eigenvalues cja > 0. We organize them as 
follows with the positive eigenvalues {uox > 0) in the first 
2N entries and the negative ones {—ujx < 0) in the second 
27V: 



cja</ca for l<iy<2N 

^-^^a$;-^^/cA for 2N^l<u<AN 





for 1< A < 2^" 
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,u-2N 



a 



-2Ar* 



for 
for 



for 27V 



-i< x<m. 



l<u <2N 



(A3) 



The eigenvectors xx form a complete and orthonormal 
set. The constants cx = cx-2N are chosen to provide the 
needed normahzation of the ax eigenvectors, which are 
not orthogonal, and are not normalized to unit norm, as 
described below. 

Because the linear eigenvalue problem is hermitian, 
its eigenvectors are orthonormal, implying x^x\, = 
5xx>' If we choose A and A' to both lie between 1 and 
2N, we find 

2N 

J2{^x^x' + ^fX^v = cxSxy. (A4) 

We choose cx = ujx/hm in order to satisfy the normaliza- 
tion condition for the collective mode raising and lower- 
ing operators [aA,a]^/] = Sxx' in Eq. (36). We will choose 



the exact same positive value for cx for the corresponding 
negative eigenvalue solutions as well {cx = cx-2n)- 



We can derive three more relations between the a 
eigenvectors by using the completeness of the x eigen- 



vectors Xlt=i^A*^A' ~ We start with v and v' 

both restricted between 1 and 2N . This gives 

^aK*c.a + = nm5,,>. (A5) 

A:u;a>0 

Next, taking I < v <2N d.nd 2N ^ I < v' < AN yields 
^ K*<'-«'*)=0. (A6) 

A:u;a>0 

Finally, taking both v and v' to be larger than 2N gives 
^ -l(ar<'+«^«A'*) = z-^<5..'. (A7) 



A:u;a>0 



^A 



In all of these relations, the second terms in the paren- 
theses arise from the negative eigenvalue solutions of the 
QEP. 
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